source("funkce.R")

dfkrych = read.csv("Data\\cohp-Bulk.csv",header=T, sep=",", dec=".")
dfstab = read.csv("Data\\cohp-GB-stable.csv",header=T, sep=",", dec=".")
dfnestab = read.csv("Data\\cohp-GB-unstable.csv",header=T, sep=",", dec=".")

skrych = read.table("soubory.dat\\GB-bulk.txt")
sstab = read.table("soubory.dat\\GB-stable.txt")
snestab = read.table("soubory.dat\\GB-unstable.txt")

mkrych = read.table("soubory.dat\\GB-bulk_mat.txt") * 3.65
mstab = read.table("soubory.dat\\GB-stable_mat.txt")
mnestab = read.table("soubory.dat\\GB-unstable_mat.txt")

rownames(skrych) = c("Si1", paste("Ni", 2:4, sep = ""))
rownames(sstab) = c(paste("Si", 1:16, sep = ""), paste("Ni", 17:64, sep = ""))
rownames(snestab) = c(paste("Si", 1:16, sep = ""), paste("Ni", 17:64, sep = ""))

# bulk --------------------------------------------------------------------

k=skrych
for (i in 1:3){
  for (j in c(1,-1)){
    temp = skrych
    temp[,i] = temp[,i]+j
    k=rbind(k, temp)
  }
}
skrych = as.matrix(k) %*% as.matrix(mkrych)

dvoj = dvojice(dfkrych)

dvojr = matrix(0, nrow = nrow(dvoj), ncol = ncol(dvoj))
for (i in 1:nrow(dvojr)){
  for (j in 1:ncol(dvojr))
    dvojr[i,j] = which(dvoj[i,j] == rownames(skrych))
}

o = NULL
for (i in 1:nrow(dvojr)){
  p = voronoi(as.matrix(skrych), dvojr[i,1], dvojr[i,2])
  o = c(o, obsah(p))
  #ocoph = obspodgr(dfkrych$Energy,dfkrych[,1+i])
}



# stab --------------------------------------------------------------------

k=sstab
for (j in c(1,-1)){
    temp = sstab
    temp[,3] = temp[,3]+j
    k=rbind(k, temp)
    temp = sstab
    temp[,1] = temp[,1]+j
    k=rbind(k, temp)
}
sstab = as.matrix(k) %*% as.matrix(mstab)

dvoj = dvojice(dfstab)

dvojr = matrix(0, nrow = nrow(dvoj), ncol = ncol(dvoj))
for (i in 1:nrow(dvojr)){
  for (j in 1:ncol(dvojr))
    dvojr[i,j] = which(dvoj[i,j] == rownames(sstab))
}

o = NULL
ocoph = NULL
for (i in 1:nrow(dvojr)){
  print(i)
  p = voronoi(as.matrix(sstab), dvojr[i,1], dvojr[i,2])
  o = rbind(o, c(i,obsah(p)))
  ocoph = rbind(ocoph,c(i,obspodgr(dfstab$Energy,dfstab[,1+i])))
  #print(paste("o = ", o[length(o)]))
}

infostab = cbind(dvoj,ovor = o[,2],ocop = ocoph[,2])
info = infostab
#info = (read.csv(file = "ulozeno\\stab.csv"))
###############################################################################

info = info[abs(info$ovor) != Inf,]

info = cbind(info, podil = info$ovor/info$ocop)

info = info[info$ovor != 0,]

#write.csv(info, file = "ulozeno\\obsahy_stab.csv")

#info = (read.csv(file = "ulozeno\\obsahy_stab.csv"))
info = info[,-1]
#info = info[,-6]

plot(info$podil, main = "podil, stab")

sisi = info[substr(info$x,1,2)=="Si" & substr(info$y,1,2)=="Si",]
nini = info[substr(info$x,1,2)=="Ni" & substr(info$y,1,2)=="Ni",]
sini = info[(substr(info$x,1,2)=="Si" & substr(info$y,1,2)=="Ni") | (substr(info$x,1,2)=="Ni" & substr(info$y,1,2)=="Si"),]

plot(sini$ovor,sini$ocohp, main = "sini stab")
plot(sisi$ovor,sisi$ocohp, main = "sisi stab")
plot(nini$ovor,nini$ocohp, main = "nini stab")

plot(sini$podil, main = "podil, sini, stab")
plot(sisi$podil, main = "podil,sisi, stab")
plot(nini$podil, main = "podil, nini, stab")

sini[sini$podil < -2,]
sisi[sisi$podil < (-100),]
sini[sini$podil < -2,]


# nestab --------------------------------------------------------------------

k=snestab
for (j in c(1,-1)){
  temp = snestab
  temp[,3] = temp[,3]+j
  k=rbind(k, temp)
  temp = sstab
  temp[,1] = temp[,1]+j
  k=rbind(k, temp)
}
snestab = as.matrix(k) %*% as.matrix(mnestab)

dvoj = dvojice(dfnestab)

dvojr = matrix(0, nrow = nrow(dvoj), ncol = ncol(dvoj))
for (i in 1:nrow(dvojr)){
  for (j in 1:ncol(dvojr))
    dvojr[i,j] = which(dvoj[i,j] == rownames(snestab))
}

o = NULL
ocoph = NULL
for (i in 1:nrow(dvojr)){
  print(i)
  p = voronoi(as.matrix(snestab), dvojr[i,1], dvojr[i,2])
  o = rbind(o, c(i,obsah(p)))
  ocoph = rbind(ocoph,c(i,obspodgr(dfnestab$Energy,dfnestab[,1+i])))
}


infonestab = cbind(dvoj,ovor = o[,2],ocop = ocoph[,2])
info = infonestab

info = (read.csv(file = "ulozeno\\nestab.csv"))
############################
###################################################

info = info[abs(info$ovor) != Inf,]

info = cbind(info, podil = info$ovor/info$ocohp)

info = info[info$ovor != 0,]

info = info[,-1]
#info = info[,-6]

plot(info$podil, main = "podil, stab")

sisi = info[substr(info$x,1,2)=="Si" & substr(info$y,1,2)=="Si",]
nini = info[substr(info$x,1,2)=="Ni" & substr(info$y,1,2)=="Ni",]
sini = info[(substr(info$x,1,2)=="Si" & substr(info$y,1,2)=="Ni") | (substr(info$x,1,2)=="Ni" & substr(info$y,1,2)=="Si"),]

sisi
sini = sini[sini$ovor<0.15,]
nini = nini[nini$ovor<0.3,]

plot(sini$ovor,sini$ocohp, main = "sini stab")
plot(sisi$ovor,sisi$ocohp, main = "sisi stab")
plot(nini$ovor,nini$ocohp, main = "nini stab")

plot(sini$podil, main = "podil, sini, stab")
plot(sisi$podil, main = "podil,sisi, stab")
plot(nini$podil, main = "podil, nini, stab")

sini[sini$podil > -0.002,]
sisi[sisi$podil < (-100),]
nini[nini$podil > (-0.05),]

sisi[sisi$ovor>20,]

